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ABSTRACT 

We present a potential-density pair designed to model nearly isothermal star clusters (and similar 
self-gravitating systems) with a central core and an outer turnover radius, beyond which density falls 
off as r~ 4 . In the intermediate zone, the profile is similar to that of an isothermal sphere (density 


pocr 2 ), somewhat less steep than the |King (1962) profile, and with the advantage that many 


dynamical quantities can be written in a simple closed form. We derive new analytic expressions for 
the cluster binding energy and velocity dispersion, and apply these to create toy models for cluster 
core collapse and evaporation. We fit our projected surface brightness profiles to observed globular 
and open clusters, and find that t he quality of the fit is generally at least as good as that for the 
surface brightness profiles of King (1962). This model can be used for convenient computation of the 
dynamics and evolution of globular and nuclear star clusters. 

Subject headings: globular clusters: general — open clusters and associations: general — galaxies: 
nuclei methods: analytical 


1. INTRODUCTION 

Although stellar systems can be modeled directly from 
observations in non-parametric ways, parametrized mod¬ 
els are often desirable as sources of intuition, theoretical 
insight, or computational convenience. Clusters of stars 
have been modeled with a wide variety of parametrized 
functions; many of t hese are sum mari zed in Chapters 
2 and 4 of Binney & Tremaine (2008). The simplest 
comm only used analytic mo dels, for example |Plummer| 
(1911) and Hernquist (1990), are two parameter fits that 
do not allow for a variable ratio between outer (or half 
mass) radius and core size. Since this ratio varies, and 
can be quite large for realistic star clusters, three param¬ 
eter models generally can offer much better fits. 

In this paper we use a three-parameter potential- 
density pair which, in projection, fits observations of 
globular and open clusters qu ite welF l Although its func¬ 
tional form is remin iscent of Deh nen (1993), the surface 
density is closer to |King| (|l962|), a profile widely used 
for fitting observed cluster surface brightnesses. When 
deprojected, this model is similar to the physically moti¬ 
vated isothermal sphere but with the additional advan¬ 
tage of analytic tractability: quantities such as density, 
potential, and binding energy can easily be written in 
closed form, but canno t for the isothermal sphere or the 
related model of King (1966). Three important features 
ar e shared w ith the King profile (and the King model 
of King| | 1966): a central, constant-density core, an exte¬ 
rior modification that produces finite mass, and an in¬ 
termediate power law zone inside the halo (or, for King, 
truncation) radius. In this intermediate zone our profile 

1 After this work was submitted for publication it was pointed 
out to us that this potentia l-density pair has bee n used for modeling 


gravitational lensing, e.g. Kassiola & Kovner (1993). Aside from 
the primarily numerical study ot Merritt & Valiuri (1996), which 
examined a triaxial analogue to this paper's potential-density pair, 
they have not been used for dynamical investigations of the equi¬ 
librium and evolution of nearly isothermal stellar systems. 


scales as p a r~ 2 and the King profile as a r^ 3 . Al¬ 
though globular and nuclear star clusters are our primary 
motivation, this potential-density pair may be useful for 
fitting galactic potentials, as between the core and halo 
radii, it possesses nearly flat rotation curves. Other pos¬ 
sible applications include the core d dark matter profiles 
found in dwarf spheroidal galaxies (Burkert|2015), or the 
“/3-mod els” used to fit X-ray o bs ervatio ns ol galaxy clus¬ 
ters (Cavaliere & Fusco-Femiano 1976). 

In 1 J 2 ] we present the basic formulae for our model: 
exactly when possible, and in approximation when nec¬ 
essary. In (j3j we consider two major applications for 
our model: as a theoretical testbed for studies of cluster 


c ore c ollapse (13.1), and as a tool for fitting observations 
(j |3.2[ ). We summarize and outline future uses for this 
model in 0 

2. POTENTIAL-DENSITY PAIR 

We introduce a stellar density profile, which in func¬ 
tional form is related to the well-known Dehnen or “ 77 ” 


profiles ( Dehnen||1993 Tremaine et al.fl.994 1: 


p[r) = 


(1 + 


r 2 /r 2 )(l 


2 Ah) ‘ 


( 1 ) 


Here p c is the central density, r c is the core radius, and 
rh is the outer halo radius (which approaches the half¬ 
mass radius when r c <C r\j). This density profile has been 
chosen to produce a flat central core contained interior to 

a roughly isothermal p a r~ 2 cluster. The halo densities 

_-3.5 


due to outer anisotropy (Spitzer & Hart 

1971 

Spitzer 

& Shapiro||1972 

1 ; we have chosen the somewhat 

steeper 


par ’ scaling for analytic tractability. 

We compare this profile to the classic isothermal 
sphere solutior0 in Fig. [l] Interior to the halo radius, 

2 I.e. the solution to ^ 3 A + WF r A = 0 with a constant 

and non-singular inner boundary conditions. 
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Fig. 1. — Comparison between the classic isothermal sphere 
(solid lines) and Eq. Id with halo radii at r h = 10 1 ' 5 (dotted 
curves; indicated by arrow) and = oo (dashed curves). For 
all curves we set G = p c = 1, and show solutions for a = 1 (black) 
and cr = 2 (green). 

rh, the density profiles are quite similar, provided we 
set the core radius r c = ri so -\/2/3, where the isothermal 
sphere core radius rj so = 3a/t/4nGp c . This ensures an 
asymptotic (r c <C r <C rh) match between p(r ) and the 
isothermal sphere. 


In projection, our model’s surface brightness profile is 


m = 


rh 


T(r h - r c ) l y/l+lp/rl y/1 + i? 2 /rg 


( 2 ) 

where I? is a projected 2D radial coordinate and T is the 
mass-to-light ratio (assumed to be constant in r). We 
have defined a central surface density E c = 7rp c r c rh/(rh+ 
r c ). For R < rh, a suitable choice of par ameters can 
make Eq. [2| very simila r to a King profile (King 1962); 
we elaborate on this in i 3.2 The potential is 


$(r) = — 


4iTp c GrZrl r 




— arctan(r/rh)-- arctan (r/r c ) 

r r 


( 3 ) 


Relatedly, the mass enclosed at a radius r is 



Fig. 2.— Top panel: the velocity dispersion a for the potential- 
density pair. Here we show four clusters with Mtot = 10®M® and 
rh = 1 pc; the red, purple, blue, and black lines show r c = 0.1 pc, 
r c = 0.05 pc, r c = 0.01 pc, and r c = 0.002 pc. Dot-dashed 
lines on the right are the virially averaged cr v , while solid lines 
are approximate closed form solutions <7 n ear(r) and o"f ar (r). 
Dotted lines represent exact numerical solutions of Eq. [9] Taken 
piecewise, our approximate solutions describe velocity dispersion 
well for r c < O.lrh, but break down severely for r c > 0.2rh- 
Bottom panel: using the same line colors, dashed lines show 
circular velocity curves, v c { TC (r) = yjGM(r)/r, and solid lines 
show escape velocities u eS c(r) = y/—2<&(r) for each cluster. 


The above formulae have been used in the p ast to stud y 
gravitational lensing (e.g. appendix A of Ehasdottir 
et al. 2007); we now derive new properties of this 
potential-density pair more relevant for star cluster dy¬ 
namics. The total potential energy of the cluster is 
W = —47 tG f Q rp(r)M(r)dr. This integral evaluates 
to 


W =- 


Wq 

rh~r c l 


r c ln4 + r h ln4-2r c ln(l+r h /r c ) (8) 


M(r) 


4nr c r h pc ^ arctan ( r / rh ) 


— r c arctan(r/r c )]. 


( 4 ) 

This expression converges as r —> oo, so the total cluster 
mass is 




2 ?r 2 rgrgp c 
rh + r c 


( 5 ) 


When r c <C rh, the half-mass radius is nearly and the 
core mass is given by 


- 2r h ln(l+r c /r h ) , 

where we have defined Wq = —GM/ ot 7r -1 /(rh — r c ). By 
the virial theorem, the total cluster binding energy E = 
IK/2. 

Computation of the velocity dispersion, a(r), is more 
challenging. The virially averaged value is a v = 
(—W/Mtot) 1 / 2 , but finding the radial dependence is more 
difficult. Assuming isotropic velocities, the exact three- 
dimensional dispersion is 


M cS M(r 


ni 


( 6 ) 


The escape velocity from any point in the cluster is 
v eS c(r) = \J— 2$(r). At the cluster center, this becomes 

u esc (0) = ^ /EMfUl yi n ( rh / rc ) _ (7) 

V 71 r h - r c 


2/ \ 3G 


M(r')p(r') 


,12 


dr' 


p(r) j r r 

12irGp c rlr^ (r 2 + r 2 )(r 2 + r/) 


(9) 


rt - r± 


rt - rt 


Y(r, r h ,r c ) 


- Y(r , r c , r c ) - Y (r, r h , r h ) + Y (r, r c , r h ) 
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where 


v , \ n arctan(r'/n) , 

y(r,ri,r 2 )=/ , /2 / 2 dr ■ ( 10 ) 

J r r 4 1 + r^/r| 

Analytic approximations for cr(r) are derived in the ap¬ 
pendix for the limit of r c <C ry. In this regime, the 
central velocity dispersion is 


r 2 = lim a 2 

r->0 


6 GM to t( 7 r 2 /8 — 1) 


nr h 
„2/ 2 , 


(ID 


l2irGp c r 2 (n 2 /8 — 1). 


Interestingly, in this limit, <r 2 /cr 2 ~ 1.011. In general, 
clusters with r c ~ ry are almost perfectly isothermal 
for r < ry/ 2 , while those with smaller core radii see a 
plateau of elevated velocity dispersion for ry < r < ry. 
In this plateau the dispersion is nearly isothermal, with 
cr 2 ss 3GM tot /(irr h ) (for ry <C r <C ry). Although the 
temperature inversion visible for ry < r < ry is not phys¬ 
ical for long-term cluster equilibria, it is characteristic of 
a post-collapse cluster core undergoing gravothermal ex¬ 
pansion. Fig. [2]shows exact numerical solutions for cr(r), 
our approximate analytic <7f ar (r) and cr near (r) derived in 
Appendix A, and other velocities of interest (u c j rc , v esc ). 

The isotropic distribution function (DF) can be cal¬ 
culated numerically using Eddington’s formula, but a 
closed form does not appear to exist. Thi s is i n con - 
trast to the somewhat similar King model ( King| 1966). 
It is possible, however, to modify many of the above for¬ 
mulae into accomodating a central black hole (BH) of 
mass to,. Leaving p(r) and I(R) unchanged, and setting 
M,(r) = M(r)+m . and $,(r) = $(r) — Gm./r, we find: 


W, = W — 


47 TGm.r 2 r 2 p c ln(r h /r c ) 


er 2 (r) = a 2 (r) + 3 Gto. ( 1 + — ) ( 1 + 


( 12 ) 

(13) 


7r ( r h + ryry + r 2 ) arctan(r/ry) r 2 arctan(r/ry) 


2 ryry(ry + ry) 


r c (rl - r 2 ) 


rh{rl - r 2 ) 


It is not obvious, but <r 2 — a 2 oc 1/r at large radii due to 


multiple high-order cancellations in Eq. 13 
Under the approximation t hat ry <C ry, we calc ulate 
the two-body relaxation time (Binney fe Tremaine|2008), 


t T (r) = 


6.5 x 1CT 2 er 3 (r) 
G 2 TO*p(r) In A 


(14) 


w ith the C oulomb logarithm A « 0.4Af tot /m t (Spitzer & 
Hart||1971 ). Specifically, the central relaxation time 


0.39 

tr(0) “ h^A 1 


r 3 M t ot \frpr h 


GM U 


TO* 


r h 


(15) 


Approximating the inner regions of the cluster as isother¬ 
mal lets us write f r (r) ~ t r ( 0 )(l+r 2 /r 2 )(l-|-r 2 /r^), which 
is valid until the turnover in the velocity dispersion at 
r « r h / 2 . 

3. APPLICATIONS 
3.1. Single Component Cluster 


If we apply the potential-density pair to a cluster com¬ 
posed of single, identical stars of mass m*, neglecting 
three-body and tidal capture processes, we can describe 
very simply the dynamics of core collapse. In particular, 
under the assumption of self-similarity, we can write 


dry ry 
dt * f r ( 0 ) 


(16) 


Self-similar core collapse is expected in s ingle-mass clus - 
ters on analytic grounds (Lynden-Bell & Egg leton .1980), 
and this is verified in numerical simulations (|Baumgardt 
et ah 20031 and references therein). In the context of our 
model, self-similarity occurs under a simplifying assump¬ 
tion about heat conduction. Specifically, we assume that 
the cluster lacks internal sources of heat (i.e. no binaries, 
no BHs, etc.) and set 


E, 


cond 


= A 


cond' 


M{r c )a\ 

2f r (0) 


d E c dry 
dry d t 


(17) 


where the binding energy of the core E c « M c a 2 / 2. 
This evolution is trivially self-similar for constant den¬ 
sity cores as long as we take the core cr to be isothermal 
(a good approximation for our model), and set A con( j to 
a constant value. This latter step must unfortunately 
be done artificially; the mild “bump” in our model’s a 
profile prohibits a self-consistent calculation of outward 
heat flux. 

We find an analytic solution in the limit of ry <C ry 
and fixed ry: 


r c (t) = \/r*c (0) - 5.1U4 cond In A—— \jGM tot r h (18) 

■Mtot 


Thus, the core-collapse timescale in this limit is 
1 


T — 

c.c 


2 A, 


cond 


0.19 


-tr, O ( 0 ) 


(19) 


_ ' r-g M tot r 2 

A mnd lnAV GM tat m* r 2 ' 


^cond 

where t r ,o( 0 ) is the initial central relaxation time. 
Isotropic Fokke r-Plan ck models take « 330 central re¬ 
laxation times ( Cohn|[l98 0 ) to reach total core collapse, 
implying A cond R 3 1.5 x 10 -3 . Direct N -body simulations 
are in reasonable agreement with this (Baumgardt et al. 
20031. 

With this calibration, we can now self-consistently 
model the evolution of our cluster. Heat flowing outward 
from the collapsing core will inflate the outer regions (as¬ 
suming an isolated cluster unaffected by external tides) 
so as to satisfy dW/dt = 0. This produces the differential 
equation 


dry 

dt 


dry W - In [1(1 + ry/ry)] W 0 
d^W + ln[I(l+r c /ry)]W 0 ' 


( 20 ) 


Eqs. (16 1 and (20) provide a simple description of a 


cluster undergoing core collapse while conserving energy. 
We plot the core radii and outer radii of core collapsing 
star clusters in Fig. [3] We see that ry evolves very little, 
justifying the approximation used in Eq. [Ts} 

Competi ng against core collaps e is the evaporation of 
the cluster (Ambartsumian |1938), due to losses from the 
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purple lines correspond to r c = 0.1,0.01, and 0.001 pc, respec¬ 
tively. If we integrate these curves with respect to radius, we find 
total evaporation rates of N e sc = 1.3 X 10 —6 yr —1 ,2.2 x 10 —6 yr —1 , 
and 2.6 X 10 —6 yr _1 , respectively. 


Fig. 3.— The core radius r c (panel A) and halo radius 
(panel B) of an equal-mass star cluster undergoing core collapse. 
The self-similar process runs away rapidly, and is highly sensitive 
to initial conditions. The orange, black, and purple curves 
correspond to Mtot of 1O 4 Mq, 1O 5 M0, and 10 q Mq, respectively, 
and in each of these cases the initial conditions are r h = 1 pc, 
r c = 0.1 pc. The grey curves are the same as the orange ones, but 
with r c = 0.01 pc initially. The thic k d otte d and dashed curves are 
numerical integrations of Eqs. |17| and |20| with fixed and variable 
central cr, respectively, while the thin soli d curve represents the 
self-similar analytic solution from Eq. |18| All solutions come to 
overlap as r c < 


high-velocity tail of the velocity distribution. We assume 
a Maxwellian velocity distribution f(v) oc exp(— f^), 
and use this to define the fraction of stars in the distri¬ 
bution that will escape the cluster, 


Cesc ir ) 


r , x f(v)v 2 dv 

JVesc(r) J V > 

r f(v)v 2 dv 


fd V esc (r) - 
V 7r er(r) 



( 21 ) 


The approximate equality assumes v esc (r) > < 7 ( 7 -). As¬ 
suming that this fraction is replenished on a local re¬ 
laxation time, t T (r ), the total rate of stellar evaporation 

from the cluster is N esc = f (0 °° d ^ ac dr, where 


d N PJ 


dr 


: dr = 47T?’ 2 


Pir) Cesc(r) 
m* t T (r) 


( 22 ) 


We plot d-/V esc /dlnr in Fig. [7j Most of the mass flux 
out of the cluster comes from radii r c < r < r^. For 
fixed M to t and rh, decreasing the core radius r c will 
marginally increase the total evaporation rate N esc . Un¬ 
surprisingly, the specific evaporation rate is much higher 
outside the core than in it. If we approximate N esc ss 
dV/dlnr| r=rh , we can compute an evaporation timescale 

7 ev;i p — V/TVegc US 


, _ 14 I r h M t°t 

evap ~ In A y GMtot m* 


(23) 


in the limit where r c rh- Comparing to Eq. |19[ 

we see that for a cluster without internal heat sources, 
Tevap/T cc ~ 0.1(rh/r c ) 2 : core collapse will proceed rela¬ 
tively unimpeded by evaporation if r c < O.lrh, but sig¬ 
nificantly greater initial core radii will allow evaporation 
to dominate at early times. 


3.2. Parametrized Fitting to Observations 

Despite its pleasant analytical features, this model is 
primarily useful insofar as it can describe realistic star 
clusters. In this subsection, we compare in greater detail 
to th e well-known King surface brightness profile (King 
1962), 


Ik(R) = 


1 


1 


1 + R 2 /r 2 c K 


1 


' t.K 


Ac 


(24) 

The three free parameters of the King profile are a char¬ 
acteristic surface density E c k, the core radius r Cj k, and 
the truncation radius rt,K- The numerically computed 
King density profile differs from Eq. [I] outside a shared 
constant-density core, as they p ossess different power law 
indices (r -2 for Eq. |T| 7 ,_3 for King|l 962|. Furthermore, 
our p(r) declines oc r~ 4 beyond the halo radius, while 
the King profile falls to p = 0 at r = r t , k- Eq. Tl is 
likely more appropriate for clusters which are not tiaally 
truncated, as was seen in the application of untruncated 
power la w mo d els t o many LMC young massive clusters 
(Elson et al. 1987). Applied to galactic globular clus- 
ters, these untruncated models generally outperform the 
King model when clusters can be resolved out to many 
half-mass radii, although the advant age i s less clea r than 
it is f or young massive clusters (McLaughlin & van der| 
Marel 2005]) , and the globulars can also be we ll-tit by 
the truncated (though more extended) model of Wilson 
(19751. 


We can relate Eq. [Tl and the King profile by equating 
their central surface density, core radius, and enclosed 
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(C) 



IT 


Fig. 5.— Panels A and B show best fit parametrizations for 
both our model I(R) (brown, solid) and the King profile Ik(R) 
(green, dashed). Panel A shows the LMC globular cluster NGC 
2019, where both parametrizations fit the data well. Panel B 
shows the LMC globular cluster NGC 1916, where our model fits 
signif icantly better than a King profi l e. In b oth cases data is taken 
from |McLaughlin & van der Marel| (2005]), and one-cr error bars 
are shown (a nandlul ot data points without error bars are shown 
for completeness, but were not used in the fit). Panel C shows the 
distribution of r 2 goodness-of-fit statistics for our model (brown) 
and the King profile (green). The results are very similar. 


mass. Respectively, 
7T r h r c p c 


Th 


= £ r 


r c = r c> K 


1 


r h = 


1 - 


yfX 


y/l + X 


arctan(VX) 


2tt L1 + X 
21n(v / V + VTTV) 


v / TTx 


X 1 


(25) 

(26) 
(27) 


VT+x 


In all three equations we have used X = 1kAc 2 ,K- To¬ 
gether, the above equations translate the parameters, 


{Sc, k, r c ,K, r t,K} of a given King profile into those for 
our model, {p c ,r c ,rh}- The reverse procedure must be 
done numerically because Eq. 27 is nontrivially tran¬ 
scendental. 

In Fig. [5] we show best-fit parametrizations for both 
I(R) and I k (R). We use the sample o f 125 Milky 
Way globular clusters from Trager et al. (19951, and 
an a dditi onal 68 LMC, SMC, and Fornax clusters from 


McLaughlin & van der Marel (2005 using their V-band 
profiles). We lit every cluster’s tabulated surface bright¬ 
ness profile, weighting our fits appropriately (using the 
measurement errors at each isophote for |McLaughlin fc| 
|van der 'Marel 20051 and more ad hoc data weights pro¬ 


vided in ITrager et al. 1995). We find that within the 
sample, mean and median r 2 values for nonlinear model 
fitting are extremely similar for both I(R) and Ik(R)- 
In particular, our mean r 2 statistic for the I(R) (Ik(R)) 
profile is 0.952 (0.942), and the corresponding median 
value is 0.974 (0.969). A histogram of the r 2 distribution 
is shown in Fig. [5] 


4. CONCLUSIONS 

We have presented a novel potential-density pair that is 
cored, truncated, and approximately isothermal (within 
the halo radius). Many important properties can be writ¬ 
ten in closed form. In this paper, we have calculated the 
binding energy, surface brightness profile, central veloc¬ 
ity dispersion, central escape velocity, and central relax¬ 
ation time. The primary quantity of interest that does 
not appear reducible to closed form is the stellar DF. 

The properties of this model match, surprisingly well, 
those of observed star clusters. We have fit our surface 
brightness profile to a large number of galactic globular 
clusters and found generally high goodnesses of fit, com¬ 
parable to or slightly better than best-fit King profiles. 

Because of this model’s simplicity and analytic 
tractability, it is well-suited for the idealized study of 
cluster dynamics. In this paper we have presented a 
closed-form description of the evaporation and core col¬ 
lapse of single-component, single-particle clusters, and 
we hope to apply this model to other dynamical pro¬ 
cesses in the future. 

Finally, we note that the cluster model introduced here 
can facilitate inexpensive computation of the dynam¬ 
ics of mass segregated subclusters. In particular, this 
potential-density pair can provide a dynamic background 
potential in which to run N-body simulations of the cen¬ 
tral subcluster. We are currently using this technique to 
investigate intermediate mass black hole forma tion due 
to collisional runaways in nuclear star clusters ( Kiipper| 


et al. in prep, 2015), but it could also be applied to 


open or globular clusters. 
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APPENDIX 

VELOCITY DISPERSION 


Here we derive the approximate formulas for the velocity dispersion near (cr nea r; r <C r’h) and far (crf ar ; r r c ) from 
the center of the cluster. The integrals Y(r, r c ,r c ) and Y(r, Th,rh) from Eq. (19J) evaluate to elementary functions: 


Y(r,ri,n) = 


-1 

24 r 3 r 


47 t r? +4rr? — 127rr 2 ri + 3ir 2 r 3 + (24r 2 ri — 8 r?) arctan — — 16r 3 In (1 4 —i') — 12r 3 arctan 2 — 

\ r \ r 2 J r i 


(Al) 

Y{r, rh,r c ) and Y (r, r c , r'h) are dilogarithms, but we can find approximate solutions by Taylor expanding the arctan 
terms in the integrands of Y(r, rh,r c ) and Y (r, r c . r'h) in the limit r' <C r h. These integrands are dominated by 
contributions from r' ss r c , so this expansion is valid when r c < n,. These approximate solutions are 


!>,n,,r c ) +(_! + _! | ln 
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r 2 + ri 


V(r,r c ,r h ) w ^ 


3irr 2 r c 


7rr c 


3r 2 

2 r c arctan — ( ^5 - 1 
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(A2) 

(A3) 


Eqs. ® and ( Al|A3 l define cr near (r), the approximate a(r) when r <C and r c . We likewise derive a crf ar (r) 
valid for r r c and r c <C r h by making the approximation that r ~3> r c in Eq. ©• Specifically, 


CTfar (d = 


+ 


6 GMtotjr 2 + rg)(r 2 + r 2 ) 

7rr h 

arctan 21n(l+r 2 /r 2 ) 


7T 1 
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7T 


71 

8 6 r 3 rh r 2 r? 2 rr? 


arctan- ^ arctan r ^ 


9r 4 

z 'h 


3r 2 rh 


(A4) 


The limit is nontrivial, but cancellations of first through fourth order terms give cr 2 ar ~ 3GM to t/(5r) when r rh- 

When r c <C rt, cr nea , T (r) and df ar (r) can be joined piecewise to approximate <j(r) for all r. 



















